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Abstract. In this article we present a new family of high order accurate Arbitrary 
Lagrangian-Eulerian one-step WENO finite volume schemes for the solution of stiff 
hyperbolic balance laws. High order accuracy in space is obtained with a standard 
WENO reconstruction algorithm and high order in time is obtained using the local 
space- time discontinuous Galerkin method recently proposed in [20] . In the Lagrangian 
framework considered here, the local space-time DG predictor is based on a weak 
formulation of the governing PDE on a moving space-time element. For the space- 
time basis and test functions we use Lagrange interpolation polynomials defined by 
tensor-product Gauss-Legendre quadrature points. The moving space-time elements 
are mapped to a reference element using an isoparametric approach, i.e. the space- 
time mapping is defined by the same basis functions as the weak solution of the PDE. 
We show some computational examples in one space-dimension for non-stiff and for 
stiff balance laws, in particular for the Euler equations of compressible gas dynamics, 
for the resistive relativistic MHD equations, and for the relativistic radiation hydrody- 
namics equations. Numerical convergence results are presented for the stiff case up to 
sixth order of accuracy in space and time and for the non-stiff case up to eighth order 
of accuracy in space and time. 

Key words: Arbitrary Lagrangian-Eulerian, finite volume scheme, moving mesh, high order 
WENO reconstruction, local space-time DG predictor, moving isoparametric space-time elements, 
stiff relaxation source terms, Euler equations, resistive relativistic MHD equations, relativistic ra- 
diation hydrodynamics 



1 Introduction 



We present a new class of high order one-step Arbitrary Lagrangian-Eulerian (ALE) fi- 
nite volume schemes for stiff hyperbolic balance laws. While the mesh is fixed in an 
Eulerian description, in Lagrangian type schemes the computational mesh moves with 
the local fluid velocity. That means that material interfaces are moving together with the 
mesh and thus one can precisely identify their location. In the recent past, a lot of work 
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has been carried out to develop Lagrangian methods. Some algorithms are developed 
starting directly from the the conservative quantities such as mass, momentum and total 
energy [43, 53] while another class starts from the nonconservative form of the governing 
equations [4,6,64]. In any discrete scheme one has to decide where to place the degrees 
of freedom of each physical variable. The existing Lagrangian schemes in literature can 
be generally separated into two main classes: 1) staggered mesh methods, where the ve- 
locity is defined at the cell interfaces while the other physical variables are located at the 
cell center and 2) cell-centered methods, where all variables are defined at the cell center. 

In [46] Munz presented several different Godunov-type finite volume schemes for 
Lagrangian gas dynamics and, in particular, he was the first to introduce a Roe lineariza- 
tion for Lagrangian gas dynamics. It was found that the Lagrangian Roe linearization 
actually does not coincide with the Eulerian one [46]. The resulting maximum signal 
speeds of this Roe linearization have subsequently been used to construct robust HLL- 
type Riemann solvers in Lagrangian coordinates. Carre et al. [7] describe a cell-centered 
Godunov scheme for Lagrangian gas dynamics on general multi-dimensional unstruc- 
tured meshes. Their finite volume solver is node based and compatible with the mesh 
displacement. In [17], Despres and Mazeran propose a way of writing the equations 
of gas dynamics in Lagrangian coordinates in two dimensions as a weakly hyperbolic 
system of conservation laws. The system contains both, the physical and the geometrical 
part. Based on the symmetrization of the formulation of the physical part, the authors de- 
sign a finite volume scheme for the discretization of Lagrangian gas dynamics on moving 
meshes. In [36], Jua and Zhang present a high-order Lagrangian Runge-Kutta DG scheme 
for the discretization of two-dimensional compressible gas dynamics. The scheme uses a 
fully Lagrangian form of the gas dynamics equations and employs a new HWENO-type 
reconstruction algorithm as limiter to control the spurious oscillations, maintaining the 
compactness of RKDG methods. The time marching for the semi discrete schemes for 
the physical and geometrical variables is implemented by a classical TVD Runge-Kutta 
method. The scheme has been shown to achieve second order of accuracy, both in space 
and time. Another Lagrangian discontinuous Galerkin finite element method has been 
recently proposed in [44]. The method preserves discrete conservation in the presence of 
arbitrary mesh motion and thus obeys the Geometric Conservation Law (GCL). 

In a series of articles [41-43] Maire et al. develop a general formalism to derive first 
and second order cell-centered Lagrangian schemes in multiple space dimensions and 
also on general polygonal grids. In [41] the time derivatives of the fluxes are obtained 
through the use of a node-centered solver which can be viewed as a multi-dimensional 
extension of the Generalized Riemann problem methodology introduced by Ben-Artzi 
and Falcovitz [3], Le Floch et al. [5,30] and Titarev and Toro [56,57,60]. In their re- 
cent papers [10, 40] Cheng and Shu developed a class of cell centered Lagrangian finite 
volume schemes for solving the Euler equations which are based on high order essen- 
tially non-oscillatory (ENO) reconstruction, both with Runge-Kutta and Lax-Wendroff- 
type timestepping. To our knowledge, the Lagrangian schemes developed in [10, 40] 
are the first better than second order non-oscillatory Lagrangian finite volume schemes 
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published so far. Further work of Cheng and Shu contains the construction of symmetry- 
preserving Lagrangian schemes, see [11, 12]. 

A completely different class of fully Lagrangian methods can be found in meshless 
particle schemes such as the SPH approach [26-29, 45], which has become very popular 
to simulate fluid motion in complex deforming domains due to its algorithmic simplicity 
and high versatility and flexibility. 

Apart from real Lagrangian methods, where the mesh actually moves with the lo- 
cal fluid velocity, and Arbitrary Lagrangian Eulerian (ALE) schemes, see e.g. [34,48,53], 
where the mesh moves with an arbitrary mesh velocity that may or may not coincide 
with the real fluid velocity [35], there also exist Semi-Lagrangian schemes. This kind 
of method is used in general to solve transport equations. Here, the discrete solution 
is represented on a fixed Eulerian grid. However, the solution at a mesh point at the 
new time f" +1 is computed from the known solution at time t n by following back the La- 
grangian trajectories of the fluid to the end-point of the trajectory, which in general does 
not coincide with a grid point. The unknown solution at the end-point of the Lagrangian 
trajectory is then obtained via interpolation from the known discrete solution at time t n 
at surrounding mesh points, see [8,9, 13,39,49]. 

In this article we introduce a new and better than second order accurate Lagrangian 
one-step WENO finite volume scheme for the solution of stiff and non-stiff nonlinear 
systems of hyperbolic balance laws. The high order of accuracy in space is obtained 
using a WENO reconstruction [2,20,37] and the one-step time discretization is based on a 
high order accurate predictor, for which a local space-time discontinuous Galerkin finite 
element scheme is used [20,24,33]. 

The outline of this article is as follows: in Section 2 we describe the numerical scheme, 
while in Section 3 we show numerical results for three different sets of equations, namely 
for the compressible Euler equations, for the resistive relativistic MHD equations (which 
provides a natural benchmark of stiff problems) and for the relativistic radiation hydro- 
dynamics equations. Results for shock tube problems are shown as well as numerical 
convergence results for smooth solutions to validate that the designed order of accuracy 
of our schemes is reached. Finally, in Section 4 we summarize our results and give some 
concluding remarks and an outlook concerning possible future extensions of our method. 



2 Numerical method 

In this article we consider general nonlinear systems of hyperbolic balance laws of the 
form 

^ + ^^ = S(Q), ieO(f)cR, teB+, QGQ Q CR V , (2.1) 

where Q is the vector of conserved variables, Qq is the space of admissible states (state- 
space), Q(f) is the variable spatial computational domain, F(Q) is the flux vector and 
S(Q) is a nonlinear algebraic source term, which can also be stiff. 
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2.1 ALE-Type One-Step Finite Volume Schemes 

The computational domain Q is discretized by a set of moving mesh points x i+ i =x i+ i (t) 
that move with a general mesh velocity V j+ i = V i+ i (Q i+ i/X i+ i,t), i.e. 

j- t x iH = V iH (Q i+l ,x i+l ,t). (2.2) 
The spatial control volumes are defined at the current time t" as Tf = [x n 1 ;x" 1 ], where 

i-j 1+2 

we have used the notation x" +1 = x- + i (t n ). By integration over the moving space-time 

control volume (£)] x [t n ;t n+1 ) and application of Gauss' theorem, the follow- 

ing integral formulation for the balance law (2.1) is obtained: 

Axf+ 1 Qf+ 1 = A<Qf-AtfF]l 1 -F^ 1 )+A<AfS i/ (2.3) 



with the mesh spacing at time t n Ax" = x" + 1 — x"_ x and the time step At = t" + — t n . The 
cell average at time f is defined as 



X 



Qf = i / Q(x,t n )dx, (2.4) 



Axf 



the source term is given by 



I 



tu x._i(f) 

2 



and the flux at the cell interface on the moving mesh is defined as 

r a7 / {^ + ^)'i)-^ + \m{x i+h {t) r t)) dt. (2.6) 



By choosing V i+ i =0, Eq. (2.3) reduces to a classical Eulerian finite volume scheme on a 
fixed mesh, while a Lagrangian-type finite volume scheme is obtained by choosing V i+ i to 
be the local fluid velocity. Obviously, any other choice of V i+ i is also possible and within 
the Arbitrary-Lagrangian-Eulerian (ALE) framework. For convenience of notation, we 
also introduce 

F y (Q,V)=F(Q)-yQ, and A V (Q,V) = — . (2.7) 
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While eqn. (2.3) is an exact relation, a numerical scheme is obtained by using a numerical 
flux ¥Y(q7,qf) , which is a function of two arguments, namely the states q^ = q;,(x7 1 ,t) 

and qf = q^(x + j ,t) on the left and on the right of the interface, respectively: 

1+2 

f n+l 

f ^ := a7 I *!{<i^/)Mx^/))dt. (2.8) 

In this article, we use two different numerical fluxes, either a simple Rusanov-type flux 
[52], or an Osher-type flux as introduced in [23]. The Rusanov-type flux reads 

Fr(q^q,T) = ^(F y (q,7^-^)+F y (q^^))-^max(q ? r-q / 7), (2.9) 



where s max =max(max(|A(A y (q ?! ,V i+ i )) |),max(|A(A y (q^,y ! - + i )) |) is the maximum sig- 
nal speed. The Osher-type flux according to [23] reads 




(2.10) 



where 

Y(s) = Y(<£,q+ s) = q h +s (q+ - q ;7 ) (2.11) 

is a straight-line segment path connecting the two states q^ and q£, respectively, and 
the integral in Eqn. (2.10) is evaluated numerically using appropriate high order Gauss- 
Legendre quadrature formulae (see [23] for details). In (2.10), the usual definition for the 
absolute value of a matrix holds: 

lA^RlAlR" 1 , (2.12) 

where R is the matrix of right-eigenvectors, R _1 is its inverse and |A| is the diagonal 
matrix of the absolute values of the eigenvalues of A. 

For the mesh velocity, needed in (2.2) and in the fluxes (2.9) and (2.10) we use the 
Roe-averaged velocity for Lagrangian gas dynamics according to Munz [46]: 

V i+l = l(V(q^) + V(qi)). (2.13) 

Note that in Lagrangian gas dynamics, the Roe average for the velocity is just simply 
given by the arithmetic average (see [46] for details on the derivation) and not by the 
common expression valid in Eulerian coordinates [50]. 

Finally, the new position of the mesh point x-.i at time t n+1 becomes with (2.2) and 
(2.13) 

x«+l = x1 +l+ J V i+l dt. (2.14) 
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Furthermore, in (2.3) also a discrete form of the source term S ; must be chosen. Since 
equation (2.3) only gives an evolution equation for the cell averages Q" but the interface 
flux F y j needs values at the element interface, a spatial reconstruction operator is needed 

1+2 

that produces suitable interface values from the given cell averages. The original first 
order Godunov finite volume scheme uses the simple reconstruction 

q h (x7 +1 _,t) = Q?, and q h (x+ +1 _,t) = Q? +1 . (2.15) 

Higher order spatial and temporal accuracy can be obtained by using a more sophisti- 
cated reconstruction operator, described in the following section. 

2.2 Polynomial WENO Reconstruction on Irregular Meshes 

In this paper we use the polynomial WENO reconstruction algorithm proposed in [20-22] 
that produces as output entire reconstruction polynomials and not point values at the cell 
interfaces, as the original optimal WENO scheme of Jiang and Shu [37]. Since the details 
can be found in the above-mentioned references, here we only give a brief summary 
of the algorithm supposing componentwise reconstruction in conservative variables. For 
more details on reconstruction in characteristic variables see [32,37]. The reconstruction 
polynomial of degree M is obtained componentwise by requiring integral conservation on 
a stencil 

S?= l \jT? (2.16) 

j=i-l 

with spatial extension / and r to the left and right, respectively. For odd order schemes 
there is only one central stencil (s = l), with l=r=M/2. For even order schemes, there are 
two central stencils with / =floor(M/2) + 1 and r =floor(M/2) for the first central stencil 
(s = 0) and / =floor(M/2), r =floor(M/2) +1 for the second one (s = l). For all schemes, 
the fully left-sided stencil (s = 2) has Z = M and r = and the fully right-sided stencil has 
/ = and r = M. The reconstruction polynomial for each candidate stencil is written in 
terms of some spatial basis functions tp m (0 as 

M 

w£(x,f") = £ rPmtfWm :=rp m (0^ s , (2-17) 
m=0 

with the mapping to the reference coordinate given by 

x=x(Z,i)=xli+Ax?Z, and % = £ (x r i) = ■ ( 2 - 18 ) 

Throughout this paper we use the Einstein summation convention, implying summation 
over indices appearing twice. We furthermore use the Legendre polynomials rescaled to 
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the unit interval I = [0;1] as basis functions l/? m (£). Integral conservation on all elements 
of the stencil then yields 



1+i 



l - f tp m (x)1tfdx = Qj, VTfeSf. (2.19) 



Ax; 

J v n 



and, with the definition of the primitive functions of i/> m (£) 

S 

= J rp m (Q<%, (2.20) 
o 

one obtains the following compact expression for the linear algebraic system that has to 
be solved for the unknown coefficients uflf' 



**(«*^,0)-*m(£(^ Vi?eS?. (2.21) 

Adopting the usual definitions of the oscillation indicators cr s [37] and the oscillation 
indicator matrix L/ m [20] 

W„, Mm-l^j— 3^ 3^7— -<«,/ ( 2 - 22 ) 



«=1 







the nonlinear weights o; s are defined by 

w s = 7 — — vr/ w s = =— , (2.23) 

where we use e = 10~ 14 , r = 8, A s = 1 for the one-sided stencils and A = 10 5 for the central 
stencils, according to [20,21]. The final nonlinear WENO reconstruction polynomial and 
its coefficients are then given by 

w ; ,(x,£") = t/; m (£)<, with < = X> s tv£ s . (2.24) 



2.3 Local Space-Time DG Predictor on Moving Meshes 

The reconstruction polynomials w;, (x,t n ) are then evolved locally within each element in 
order to obtain high order time accuracy. Instead of the Cauchy-Kovalewski procedure 
based on Taylor series and repeated differentiation of the governing PDE used in the 
original ENO method of Harten et al. [32], in the ADER schemes of Titarev and Toro 
[56-58,60,61] and in the Lagrangian ENO finite volume scheme with Lax-Wendroff time 



8 



discretization presented by Liu et al. [40], we use a weak formulation of the governing PDE 
in space-time based on the local space-time discontinuous Galerkin method introduced 
in [20, 24, 33], which is also capable of dealing with stiff algebraic source terms. Due 
to the element-local formulation, the method proposed here is different from the global 
space-time DG method of Van der Vegt and Van der Ven [62,63]. 

In order to get an element-local weak formulation of the PDE on the moving space- 
time control volume [x ; _i (t);x i+ i (t)] x [t n ;t n+ ], the governing PDE (2.1) is first trans- 
formed to the reference space-time element Te = [0;1] 2 . Therefore, we map the physical 
variables x and t onto the reference variables £ and x, using an isoparametric mapping, 
i.e. for the mapping of the coordinates we use the same basis functions 9 m that are also 
used to represent the numerical solution. In this article we use for 6 m the Lagrange inter- 
polation polynomials of degree M that pass through the tensor-product Gauss-Legendre 
quadrature points on the reference element Te = [0;1] 2 . For details on multidimensional 
quadrature formulae see [54]. In the following, the underlying one-dimensional Gauss- 
Legendre quadrature points and weights on the unit interval [0;1] are denoted by Q and 
ocj, respectively. Using the space-time basis functions 6\, the mapping of x and t onto £ 
and t simply reads 

x{^r) = x m 6 m ^,r), t^,r)=i m 9 m ^,r). (2.25) 

Here, m = 9 m (^,T) and the coefficients x m and i m denote the nodal coordinates in physical 
space and time and £ and T are the reference coordinates. A sketch of this isoparametric 
mapping is depicted in Fig. 1. Since the time coordinates are the same for each spatial 
node, one gets the following simple mapping for the time coordinate: 

£ = £" + A£t, (2.26) 

which reduces the Jacobian of the space-time mapping (£,t) — > (x,t) given by (2.25) to 

'=(2 £) = (»£) 

and its inverse is given by 

Now, the derivatives of the PDE (2.1) are transformed to derivatives with respect to 
the reference element using the chain rule, i.e. we get 

and with the inverse of the Jacobian of the mapping (2.28) one obtains the PDE rewritten 
in reference coordinates: 

3Q Af3F i T 5Q , . , 
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Figure 1: Sketch of a third order isoparametric space-time element. Left: physical space-time element. Right: 
reference space-time element. The interpolation nodes for the numerical solution and for the mapping, given 
by the tensor-product Gauss-Legendre quadrature points, are numbered from 1 to 9. The initial location for 
the spatial Gauss-Legendre nodes Xo,»; is a ' so highlighted. 

To simplify the notation, we introduce the following operators on the space-time refer- 
ence element Te'. 

1 11 
\f>g] r = J /(&T)g(£,TK, and (f,g) =JJ /(£,T)g(£,T)^T. (2.31) 
0+0 

For isoparametric elements, the discrete solution of PDE (2.30) is approximated using the 
same space-time basis functions 6 m that have also been used for the mapping (2.25), i.e. 

qh = qn(x / t) = O m {x,t)q m . (2.32) 

Eqn. (2.30) is now multiplied with space-time test functions 8^ (the same as the basis 
functions for the discrete solution and the mapping), and is then integrated over the 
space-time reference element Tg = [0;1] 2 . One obtains 

(i^q^ + ^q^w^ + ^^F^q^-^^At^S^qft)). (2.33) 

Here, the initial condition given by the WENO reconstruction polynomials Wj l (x,t n ) at 
time t n has been introduced in a weak form via the jump term [f^q/,— w/,] . In other 
words, the first term of (2.33) is the integral over the smooth part of the solution while 
the second term takes into account the jump of the discrete solution from t n to t n,+ . Since 
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we use a nodal basis, the interpolation polynomials for the flux and source terms are given 

h y 

F h = ¥ h (q h )=6 m (£,T)P m , and S h = S h (q h ) = 9 m (g,T)S mr (2.34) 

with 

F m = F(q m ) and S m = S(q m ). (2.35) 

Using the definitions for the WENO reconstruction polynomial (2.24) and the defini- 
tions for the discrete space-time solutions (2.32), (2.34) and (2.35) we obtain the following 
element-local nonlinear algebraic equation system: 

^Lq-+^F--4U- = f L< + A^ m S m , (2.36) 
with the definitions of the following matrices: 

'e k/ ^e m \+[9 k/ e m f\, (2.37) 



K km : 



4m = \ Q *>Tfw)' *£• = \° k '%W /' (2 ' 38) 

and 

$m=[0k>1>m] , M km = (9 k/ 9 n ). (2.39) 

Due to the nodal approach on the tensor-product Gauss-Legendre nodes, the f and t 
directions decouple and the above matrices can be easily evaluated dimension by dimen- 
sion using one-dimensional Gaussian quadrature. The system (2.36) is solved using an 
iterative method similar to the one proposed in [24,33]: 

+4A - 4m tim = Fkm^m + ^M km (2.40) 

where the stiff algebraic source term is taken implicitly (see [24] for details). A partic- 
ularly efficient strategy for obtaining an initial guess for c^ m can be found in [33]. The 
equation that determines the location of the spatial coordinates x m of the space-time ele- 
ment is 

dx 

Tt =V { x,t), (2.41) 

where V(x,t) is the local mesh velocity. In the fully Eulerian case one has V — 0, while in 
the pure Lagrangian case, V(x,t) is the local fluid velocity. For the local mesh velocity we 
use the nodal ansatz 

V h = V h (x,t) = 9 m (x,t)v m , with v m = V(x m ,t m ). (2.42) 

The initial distribution of the spatial Gauss-Legendre quadrature points at time f is 
given by 

Xo^x^+AxfU, (2-43) 
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and the spatial Lagrange interpolation polynomials passing through these points are de- 
noted by (p m . Formally a discrete version of the ODE (2.41) can then be obtained using 
again the local space-time DG method, see [18]: 

((^^7 + ^°)^ +1 = [e k ,(pm] x , m + At(e k ,9 m )v l m . (2.44) 

The weak formulation for the spatial coordinates (2.44) is iterated together with the weak 
formulation for the solution (2.40) until convergence is reached. The temporal coordi- 
nates i m are always fixed and are given by the Gauss-Legendre points and the relation 
(2.26). The space-time polynomials q/j(x,f) are computed for each element in the compu- 
tational domain and are then used as arguments for the numerical flux in Eqn. (2.8) and 
the numerical source term S, is computed using q/,(x,£) as follows: 

St = Ax^AtI I S Mx r t))dxdt. (2.45) 

' f " X. ! (t) 

' 1 

This completes the description of the high order Lagrangian one-step finite volume al- 
gorithm (2.3). 



3 Test problems 

In this section we show some computational test problems to illustrate the performance 
of the scheme in the case of the compressible Euler equations (non-stiff), the resistive rel- 
ativistic MHD equations (stiff) and for the relativistic radiation hydrodynamics equations 
(moderately stiff). 




3.1 Compressible Euler equations 

The Euler equations of compressible gas dynamics read 

pu \ 

2 + p =S(x,t), (3.1) 
u{pE + p) J 

with the fluid density p, the velocity u, the total energy density pE, a vector of source 
terms S and the fluid pressure p, given in terms of the conserved quantities by the equa- 
tion of state of an ideal gas as 

p = ( 1 -l)(pE-±pu 2 ), (3.2) 

with the ratio of specific heats 7. In this section, we define the local mesh velocity as the 
local fluid velocity, i.e. we choose 

V = u. (3.3) 
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3.1.1 Numerical convergence results 

In order to assess the accuracy of the method presented in section 2 we carry out several 
simulations of a test problem with smooth solution and exact solution on a series of suc- 
cessively refined meshes. For this purpose, we use the so-called manufactured solution 
method, which means prescribing the exact solution a priori, inserting the solution in 
PDE (3.1) and putting all terms that do not cancel into the source term S(x,t). For our 
particular test problem we choose the primitive variables of the exact solution as 

1 

p e (x,t) = l + -sin(27rx)cos(27r£), u e (x,t) =sm(2nx) cos (lizt), p e (x,t) = l. (3.4) 

From there, the vector of conserved variables can be computed as Q e = {p e ,PeU e ,pef (7 — 
1) + l/2p e Ug) T . The ratio of specific heats is chosen as 7 = 1.4. Inserting (3.4) into (3.1) 
yields the source term S(x,t). The initial computational domain is O = [— 2;2] and is dis- 
cretized by an initially uniform mesh of Nc control volumes. The boundary conditions 
are periodic. Simulations are carried out with third to eighth order Lagrangian one-step 
WENO finite volume schemes using the Osher-type flux (2.10) for one period until the fi- 
nal time t — 1.0. The Courant number is set to CFL=0.9. For each mesh the corresponding 
error in L2 norm is computed as 



and the resulting numerical convergence rates are listed in Table 1. From the presented 
results we can conclude that the designed order of accuracy of the scheme is reached. 

3.1.2 Shock tube problems 

In this section we solve a set of several shock-tube problems given in [59] and [10,40]. 
The initial conditions of the Riemann problems are 



where the initial states left and right are summarized in Table 2. The initial computa- 
tional domain is O = [xj^Xr] and is discretized with 200 equidistant cells, apart from RP5 
(Leblanc shock tube), for which 2000 cells are used, and RP0, for which only 100 cells are 
used. Simulations have been carried out with the fifth order version of our Lagrangian 
one-step WENO finite volume schemes. The first problem (RP0) is the advection of an 
isolated moving contact wave with constant pressure and velocity. Any Riemann solver 
that resolves exactly stationary contact waves in the Eulerian case should preserve ex- 
actly isolated moving contact waves in the Lagrangian case. From Fig. 2 it becomes evi- 
dent that the Osher-type flux (2.10) solves the problem exactly, without any intermediate 




(3.5) 




(3.6) 
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Table 1: Numerical convergence results for the compressible Euler equations using the third to eighth order 
version of the Lagrangian one-step WENO finite volume schemes presented in this article. The error norms 
refer to the variable p (density) at time f = 1.0. 
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points in the contact wave, whereas the Rusanov flux (2.9) adds significant numerical dif- 
fusion to the problem, as expected. For this problem, the Osher-type flux leads to a pure 
Lagrangian scheme, where the mass in each moving control volume remains constant. 

In Figs. 3-5 we show the exact solution together with the computational results 
obtained with Osher-type flux (2.10) and the Rusanov-type flux (2.9) for the other shock 
tube problems RP1-RP5. Overall, a very good agreement is noted between the numerical 
solution and the exact solution. The Osher flux resolves the contact wave very well in 
general. However, some intermediate points are generated, since in the initial phase 
of the Riemann problem waves of different nature (shock and rarefaction) overlap with 
the contact wave, thus leading to some amount of numerical diffusion in the contact 
wave. These results are as expected, since in the present paper an ALE-type approach 
has been used, which does in general not impose constant mass in each control volume, 
in contrast to purely Lagrangian schemes as the one presented, for example, in [7]. For 
the Leblanc problem, it can be easily noted that the Rusanov flux is much more robust for 
this problem due to its larger numerical diffusion compared to the Osher-type scheme. 
However, the results presented here are similar to the ones presented in [10,40]. 

3.2 Resistive relativistic MHD equations (RRMHD) 

The resistive relativistic MHD (RRMHD) equations constitute a hyperbolic system of bal- 
ance laws which has a source term that may become stiff. The equations include five 
equations for the fluid motion (conservation of mass, momentum and energy), plus six 
equations for the evolution of the electric and of the magnetic field (Maxwell equations). 
Furthermore, two additional equations are needed to maintain the constraints on the di- 
vergence of the electric and of the magnetic field. In this paper, we use the hyperbolic 
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Table 2: Initial states left and right for density p, velocity u and pressure p for the Riemann problems solved 
for the compressible Euler equations. The initial position of the discontinuity (x ( ;) and the initial computational 
domain n=[z^;XK] are also specified. In all cases 7 = 1.4, apart for RP5, where 7 = 5/3. 



RP 


PL 


ml 


PL 


PR 


u R 


PR 


Xd 




XR 





1.0 


1.0 


1.0 


0.1 


1.0 


1.0 


0.0 


-0.5 


0.5 


1 


1.0 


0.0 


1.0 


0.125 


0.0 


0.1 


0.0 


-1.0 


1.0 


2 


0.445 


0.698 


3.528 


0.5 


0.0 


0.571 


0.0 


-0.5 


0.5 


3 


1.0 


0.0 


1000 


1.0 


0.0 


0.01 


0.1 


-0.5 


0.5 


4 


5.99924 


19.5975 


460.894 


5.99242 


-6.19633 


46.095 


0.0 


-1.0 


1.0 


5 


1.0 


0.0 


0.1(7-1) 


10~ 3 


0.0 


10- 10 (7-1) 


3.0 


0.0 


9.0 



Exact solution 
WENO-FV 05 (Osher) 
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0.8 - 
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0.3 - 
0.2 - 
0.1 - 
Of" 



Exact solution 
WENO-FV 05 (Rusanov) 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 



Figure 2: Exact and numerical solution obtained with third order Lagrangian one-step WENO finite volume 
schemes for RP0 (isolated moving contact wave) at f = 0.5 using 100 cells. Left: Osher-type flux (2.10). Right: 
Rusanov-type flux (2.9). 



divergence cleaning approach according to Dedner et al. [14]. The last equation expresses 
the conservation of the total charge. In Cartesian coordinates, using the abbreviations 
9 t = Jj and 9, = the resistive relativistic MHD equations can be written as follows [24]: 
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Figure 3: Exact and numerical solution obtained with fifth order Lagrangian one-step WENO finite volume 
schemes for RP1 (Sod problem) at t = 0A (top) and RP2 (Lax problem) at f = 0.1 (bottom). Left: Osher-type 
flux (2.10). Right: Rusanov-type flux (2.9). 
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-0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Figure 4: Exact and numerical solution obtained with fifth order Lagrangian one-step WENO finite volume 
schemes for RP3 at t = 0.012 (top) and RP4 at t = 0.035 (bottom). Left: Osher-type flux (2.10). Right: 
Rusanov-type flux (2.9). 
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X 



Figure 5: Exact and numerical solution obtained with fifth order Lagrangian one-step WENO finite volume 
schemes for RP5 (Leblanc problem) at f = 6.0. Left: Osher-type flux (2.10). Right: Rusanov-type flux (2.9). 



d t D+d i (Dv i ) = 0, (3.7) 

d t Sj+diZi = 0, (3.8) 

dtr+djS^O, (3.9) 

d t E i -e'' k d j B k +d i Y = -f, (3.10) 

d^+e^djEk+d^ = 0, (3.11) 

dtY+diE^pc-KY, (3.12) 

3^+90 = -%®, (3.13) 

dtpc+dJ^O, (3.14) 



where the conservative variables of the fluid are 



D = pW, (3.15) 
S ! ' = cvW^+e^EjBk, (3.16) 
T = cvW 2 -p+l(E 2 + B 2 ) , (3.17) 

expressing, respectively, the relativistic mass density, the momentum density and the to- 
tal energy density. The spatial tensor Zj in (3.8), representing the momentum flux density, 
is 

Z'j = cvW 2 v i v j -E'E j -B i B j + [p + \(E 2 + B 2 )] 5), (3.18) 
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where Sj is the Kronecker delta, while W = 1 / \/\ — v 1 is the Lorentz factor of the fluid. In 
this paper we have assumed the equation of state of an ideal gas, namely 

p=( 1 -l)pe = 7l (co-p), (3.19) 

where 7 is the adiabatic index, 71 = (7— e is the specific internal energy and to = 
pe+p+p is the enthalpy. The source term / appearing in (3.10) is the current vector, given 
by Ohm's law, for which we assume the following expression [38,47], 

J i =p c v i +crW[E i +(&VjB k -(Elvj)j] , (3.20) 

where p c is the charge density in the laboratory frame. 

The system of equations (3.7)-(3.14) is written as a hyperbolic system of balance laws 
as in (2.1) and it has source terms in the three equations (3.10) that are potentially stiff, 
see [47]. In the stiff limit case (cr — > 00) the resistive relativistic MHD equations reduce 
to the ideal relativistic MHD equations (RMHD), for which several test problems with 
exact solution are known, see [1,31,65]. For the system (3.7)-(3.14) a family of high order 
one-step schemes in Eulerian coordinates has been proposed in [24], while in this paper 
we use a Lagrangian method. 

3.2.1 Numerical convergence results 

The smooth unsteady test case with exact analytical solution used here was introduced 
for the ideal relativistic MHD equations by Del Zanna et al. [16] and was solved for the 
first time on unstructured triangular meshes with high order PnPm schemes in [19]. Since 
the resistive MHD equations tend asymptotically to the ideal ones in the stiff limit (cr — > 
00), this is an ideal test case to assess the accuracy of our scheme in the stiff limit of the 
governing PDE system. 

The test case consists of a periodic Alfven wave whose initial condition at t = is chosen 
to be p = p = l, B { = B {l,cos(kx),sm{kx)) T , v l = -v A /B ■ (0,B y ,B z ) T , E { = -e^VjBj, and 
<p — xp — q = 0. We furthermore use the parameters k = In, 7 = | and Bo = 1/ hence the 
advection speed of the Alfven wave in x-direction is v& = 0.38196601125, see [16] for a 
closed analytical expression for Va- The computational domain is O = [0;1] with periodic 
boundary conditions and the final time at which we compare the exact solution with the 
numerical one is chosen as t = 0.5. Since in this test case the fluid velocity in x-direction 
is v x = 0, we move the mesh artificially with the fluid velocity in y-direction, i.e. we 
set V = Vy. Since this test case was constructed for the ideal relativistic MHD equations, 
we have to use a very high value for the conductivity (a = 10 8 ) in the resistive case to 
reproduce the ideal equations asymptotically. In all our computations a constant Courant 
number of CFL = 0.5 is used. 

Table 3 shows the errors and the orders of convergence measured in the L 2 norm for the 
flow variables v y and E y . The number Ng denotes the number of grid points along the 
x-axis. We find that the nominal order of accuracy M+l has been reached for all schemes 
from third to sixth order of accuracy in space and time under consideration, even for the 
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Table 3: Numerical convergence results for the stiff limit (u=10 ) of the resistive relativistic MHD equations 
(RRMHD) using third to sixth order Lagrangian one-step WENO finite volume schemes. The error norms refer 
to the variable v y (velocity in y-direction) and to the relaxed variable £„ (electrical field in y-direction). 



N G 


L 2 


0(L 2 ) 


L 2 


0(L 2 ) 




Variable v y 




Variable E y 


03 


100 


8.2658E-02 




1.1382E-02 




200 


1.2933E-02 


2.7 


1.7733E-03 


2.7 


400 


1.6965E-03 


2.9 


2.3318E-04 


2.9 


800 


2.1272E-04 


3.0 


2.9843E-05 


3.0 


OA 


100 


1.5487E-02 




3.0653E-03 




200 


5.1306E-04 


4.9 


9.7284E-05 


5.0 


400 


1.9320E-05 


4.7 


3.2709E-06 


4.9 


800 


9.3922E-07 


4.4 


1.8739E-07 


4.1 


05 


100 


5.9185E-03 




7.2904E-04 




200 


2.9331E-04 


4.3 


3.6769E-05 


4.3 


300 


4.2396E-05 


4.8 


5.3220E-06 


4.8 


400 


1.0396E-05 


4.9 


1.3026E-06 


4.9 


06 


50 


1.8839E-02 




3.6056E-03 




100 


6.2118E-04 


4.9 


1.1061E-04 


5.0 


200 


1.0376E-05 


5.9 


1.6986E-06 


6.0 


300 


8.4098E-07 


6.2 


1.5450E-07 


5.9 



electric field E y , which is one of the variables onto which the stiff relaxation source term 
is acting. 

3.2.2 Shock tube problems 

In this section we solve two out of a series of test problems proposed by Balsara in [1] for 
the ideal relativistic MHD equations. In particular, we solve the resistive RMHD equa- 
tions with a large value for the conductivity a to validate the behaviour of our method in 
the presence of stiff source terms. The initial condition is given by two piece wise constant 
states separated by a discontinuity at x=0. The left and right values for the primitive vari- 
ables are reported in Table 4. Furthermore, we set E 1 = —e^v^By (p = ip = q = and 7 = §■ 
The conductivity in our test cases is chosen as c= 10 3 for the first test problem and cr= 10 5 
for the second one. The computational domain is Q = [— 0.5;0.5] with Dirichlet bound- 
aries consistent with the initial condition in ^-direction. We use an initially equidistant 
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Table 4: Initial left (L) and right (R) states for the resistive relativistic MHD shock tube problems and final 
times t e . 



Case 


P 


V 


u 


V 


w 


By 


B z 


B x 


t e 


1L 


1.0 


1.0 


0.0 


0.0 


0.0 


1.0 


0.0 


0.5 


0.4 


1R 


0.125 


0.1 


0.0 


0.0 


0.0 


-1.0 


0.0 


0.5 




2L 


1.08 


0.95 


0.4 


0.3 


0.2 


0.3 


0.3 


2.0 


0.55 


2R 


1.0 


1.0 


-0.45 


-0.2 


0.2 


-0.7 


0.5 


2.0 





grid with 400 points. The numerical results are shown together with the exact solution in 
Figures 6 and 7. The exact solution is the one for the ideal RMHD equations, as published 
in [31]. The essential wave structures of the ideal RMHD Riemann problem can be noted. 



3.3 Relativistic Radiation Hydrodynamics 

As an additional test, we have considered the solution of the special relativistic radiation 
hydrodynamics equations. In the truncated moment formalism introduced by [55], it is 
possible to write such equations, at least in the optically thick regime, in the conservative 
form required by Eq. (2.1) (see [25], [51] and [66]). In one spatial dimension, the vectors 
of the conservative variables in eqn. (2.1), of the fluxes and of the sources are respectively 
given by 

( vD\ 
Z 
S 

Rr 
Sr 



( D \ 



Q 



D 

S 

T 
S r 

Tr 



V Tr J 



( 



\ 





Gr 

G< 

-G r 

-a 



\ 



(3.21) 



/ 



The first three equations express the usual conservation of mass, momentum and energy 
of the fluid, and the corresponding conservative variables (D,S,t) have the same defi- 
nition as in (3.15)-(3.17), except for the fact that the electromagnetic fields are zero. The 
last two equations, on the other hand, represent the time evolution of the flux and of the 
energy density of the radiation field as measured in the laboratory frame, with 



Sr = -ErW 2 V + Wfr{lW), 
Tr = h r W 2 + 2WfrV-^-. 



(3.22) 
(3.23) 



The primitive variables of the radiation field are the flux f T and the energy density E r 
as measured in the comoving frame of the fluid, and they are formally related to the 
specific intensity of the radiation l v [66]. Finally, the quantities Z and R r entering (3.21) 




Figure 6: Numerical solution obtained for shock tube problem 1 (u = 10 ) of the resistive relativistic MHD 
equations using a third order Lagrangian one-step WENO finite volume scheme and exact solution of the ideal 
RMHD equations. Results are shown for density p, transverse velocity v, pressure p and magnetic field By. 
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Figure 7: Numerical solution obtained for shock tube problem 2 (c = 10 5 ) of the resistive relativistic MHD 
equations using a third order Lagrangian one-step WENO finite volume scheme and exact solution of the ideal 
RMHD equations. Results are shown for density p, transverse velocity v, pressure p and magnetic field By. 
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Table 5: Initial left (L) and right (R) states for the two relativistic radiation hydrodynamics shock tube problems 
considered. 



Case 


P 


V 


V 




x x lp 


7 


^rad 


1L 


1.0 


4.0xl0~ 3 


0.2425 


2.0x10^ 


0.2 


5/3 


7.812 xlO 4 


1R 


3.11 


4.512xl0" 2 


8.014xl0" 2 


3.46 xl0~ 3 








2L 


1.0 


60.0 


0.995 


2.0 


0.3 


2 


1.543 xlO" 7 


2R 


8.0 


2.34 xlO 3 


0.781 


1.14xl0 3 









are defined as 

Z = wWV + p, (3.24) 
R, = -E T W 2 v 2 +2Wf T v+-±. (3.25) 

The sources of the radiation field G r and G r depend on the physical interaction between 
radiation and matter and can be written as 

G* = x'i^-^W+ix'+xlvfr (3.26) 
G r = tf(E I -47tB)vW+tf+x!')f I/ (3.27) 

where 47rB = a ra> $T^ is the equilibrium black body intensity, while x* and X s are the ther- 
mal and the scattering opacity coefficients, respectively which are ultimately responsible 
for the stiffness of these equations. The temperature T of the fluid is computed from 
the ideal-gas equation of state through the simple relation T = p/p. We also recall that, 
while the conversion from the purely hydrodynamical conservative variables to the cor- 
responding primitive variables is not analytic and it requires the solution of an algebraic 
equation [15], the conversion from the conservative radiation variables (S r ,T r ) to the cor- 
responding primitive variables (f T ,E T ) is just linear and follows directly from (3.22)-(3.23). 

In the verification of our numerical scheme, we have considered two shock-tube tests, 
with initial conditions reported in Tab. (5). The first test involves the propagation of a 
mildly relativistic strong shock, while the second one generates a smooth highly rela- 
tivistic wave. Each test is evolved in time until stationarity is reached. The semi-analytic 
solution that is used for comparison with the numerical one has been obtained following 
the strategy described by [25]. The scattering opacity coefficient x s nas been set to zero, 
while the value of the thermal opacity coefficient x i is reported in Table 5 and it produces 
configurations that are moderately stiff. Figure 8 shows the comparison of the numerical 
solution with the exact one, where we have adopted a third-order WENO reconstruction, 
with CFL = 0.4 and 100 initially equidistant grid points. In Figure 9 we show the evolu- 
tion of the residual of the density p. It can be noted that initially the residual stagnates at 
a rather high level due to the presence of many transient waves inside the computational 
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Figure 8: Numerical solution obtained for shock tube problem 1 (left panels) and for shock tube problem 2 
(right panels) of the relativistic radiation hydrodynamics equations using a third order Lagrangian one-step 
WENO finite volume scheme. Results are shown for density p, velocity v, and energy density of the radiation 
field E r . 
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Figure 9: Evolution of the residual in the variable p for test problem 1 (left) and test problem 2 (right). 



domain. Once the transient waves have left, the residual drops very quickly to machine 
precision. These results confirm the ability of the new scheme in solving the relativistic 
radiation hydrodynamics equations. It has furthermore been shown that the proposed 
high order numerical scheme is able to simulate correctly, both, time dependent prob- 
lems as well as steady state problems. 



4 Conclusions 



We have developed a new high order one-step Arbitrary-Lagrangian-Eulerian (ALE) 
WENO finite volume scheme for the solution of nonlinear systems of hyperbolic balance 
laws with stiff source terms. The presented approach has been validated against exact 
reference solutions available for smooth and discontinuous solutions of three different 
hyperbolic systems, namely the Euler equations of gas dynamics, the resistive relativis- 
tic MHD equations and the relativistic radiation-hydrodynamics equations. In all cases 
the algorithm was found to be very robust and at the same time very accurate. To our 
knowledge, it is the first time that such high orders have been reached with ALE-type fi- 
nite volume schemes. In the near future we plan to extend the schemes presented in this 
article to structured and unstructured meshes in multiple space dimensions. Again, the 
building blocks will be a high order WENO reconstruction [21,22] and a local space-time 
DG predictor [20,24,33]. 
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